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Particle-in-cell (PIC) simulations are widely used as a tool to investigate instabilities that develop 
between a collisionless plasma and beams of charged particles. However, even on contemporary 
supercomputers, it is not always possible to resolve the ion dynamics in more than one spatial 
dimension with such simulations. The ion mass is thus reduced below 1836 electron masses, which 
can affect the plasma dynamics during the initial exponential growth phase of the instability and 
during the subsequent nonlinear saturation. The goal of this article is to assess how far the electron to 
ion mass ratio can be increased, without changing qualitatively the physics. It is first demonstrated 
that there can be no exact similarity law, which balances a change of the mass ratio with that 
of another plasma parameter, leaving the physics unchanged. Restricting then the analysis to the 
linear phase, a criterion allowing to define a maximum ratio is explicated in terms of the hierarchy 
of the linear unstable modes. The criterion is applied to the case of a relativistic electron beam 
crossing an unmagnetized electron-ion plasma. 



I. INTRODUCTION 

The dynamics of collision-less plasma far from its equi- 
librium is frequently examined with particle-in-cell (PIC) 
simulations. The unique capability of PIC codes to model 
such systems from first principles on macroscopic scales 
implies that they can bridge the gap between theory and 
experiment. For example, just a few years ago it was 
still unclear if relativistic shocks exist. It was not known 
whether the motion energy could be dissipated rapidly 
enough to sustain the shock discontinuity p], Q. Such 
shocks have not yet been observed directly, because they 
do not exist in solar system plasma. Recent PIC simula- 
tions could shed light on how they develop in response to 
the filamentation (Weibel) instability |3|— ISI] - The particle 
acceleration and the generation of electromagnetic radi- 
ation within the context of active galactic nuclei |6M9L 
sup ernova remnants [TO, [H| or gamma ray bursts jl2|- 
fl5j have also been investigated. PIC simulations are 
now instrumental in investiga ting the plasma thermal- 
ization within solar flares [16L LL7j and the dynamics of 
magnetic reconnection [l8l Il9|. On a completely differ- 
ent length and density scale, the fast ignition scenario 
for inertial confinement fusion [20| has prompted within 
the last decades many numerical works focusing on the 
propagation of charged particle beams in a collisionless 
plasma [2l|,[22|. 

Current simulations employ billions of computational 
particles, placing physically realistic PIC simulations 
within our reach 123] . The inclusion of ions in PIC sim- 
ulations nevertheless remains a formidable challenge. As 
long as the system under scrutiny involves only electrons 
and positrons with the mass m e , the time scale that must 
be resolved is typically the inverse electronic plasma fre- 
quency uj^ 1 oc ^Jm~ e . Running the simulation for hun- 
dreds or thousands uo~ l captures the evolution of the 



system way beyond its linear phase. Mobile protons or 
ions in the simulation result in an additional and much 
longer timescale. A PIC simulation must then resolve 
many inverse proton plasma frequencies uj~ x oc ^/M p 

and cover a time interval that is F p = ^/M p /m e ~ 42 
times longer, if the plasma is unmagnetized and if pro- 
tons are the only ion species. A further penalty is intro- 
duced by the larger spatial scales of the ion structures. 
The size of the ion filaments is, for example, comparable 
to the ion skin depth c/uo pi while that of the electron fil- 
aments is ~ c/uo e . In principle, the simulation box size 
that is necessary to model electron-proton plasmas in- 
creases compared to that required by leptonic plasmas 
by a factor ~ , where D is the number of resolved 
spatial dimensions. 

For this reason and until now, PIC simulations that use 
the correct electron-to-proton mass ratio are restricted 
to ID systems [§[ and to 2D simulations that resolve a 
limited spatio-temporal domain [24[, while 2D PIC sim- 
ulations that cover a large domain with regard to the 
ion scales or even 3D PIC simulations normally resort 
to reduced ion masses between 10-100 electron masses 
d, [iH, [26|. Ion masses of up to 1000 electron masses 
have been used in a 2.5D simulation [4[ thanks to a low 
number of particles per cell, which is beneficial for the 
scalability of a domain-decomposed PIC code. Multi- 
dimensional plasma simulations that employ the correct 
mass ratio and capture the largest ion scales are possi- 
ble, if the electron dynamics does not have to be resolved 
accurately. Implicit PIC schemes can dissipate away the 
energy contained in the smallest scales in a form of Lan- 
dau damping [27], . The cell size can then be increased 
beyond the plasma Debye length, without restricting the 
physical accuracy of the large-scale dynamics. However, 
if both the electron and the ion dynamics must be re- 
solved simultaneously, the implicit PIC codes are equally 
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costly as the explicit ones. The speeding up of the plasma 
evolution through a reduced ion mass will remain a ne- 
cessity in particular for 3D simulations. It is thus im- 
portant to assess how the plasma dynamics changes with 
this parameter. 

Various studies exist that demonstrate the importance 
of the electron-to-ion mass ratio for the plasma dynamics 
in several types of plasma processes. Parametric studies 
of plasma shocks have addressed this issue both in the 
non-relativistic (29l43lj and relativist ic Q regimes. Sim- 
ulation studies of the interplay between electron phase 
space holes with the ions and its dependence on the mass 
ratio can be found in [32|, [33| . The impact of the mass 
ratio on the reconnection of magnetic field lines and the 
associated particle acceleration has been investigated in 
Ref. [34],[35j. However, these studies related to the effects 
of a reduced ion mass focus primarily on the nonlinear 
evolution of the simulation. 

This article is a first systematic study of the conse- 
quences of a reduced ion mass within a theoretical frame- 
work. The impact of the mass ratio on the nonlinear cou- 
pling of the plasma dynamics across the different scales is 
not considered here; it is too complex and multifaceted. 
An example would be the enlargement of the foreshock of 
a perpendicular shock with an increasing ion mass, which 
influences the resulting instabilities and the thermaliza- 
tion of the shock-reflected ion beam [29l43ll l36| . We 
study here the spectrum of linearly unstable waves, which 
should depend on the mass ratio between the ions and the 
electrons. Ions only one time "heavier" than electrons 
are obviously too light as they behave like positrons, not 
like protons. Is it therefore possible to draw a line, from 
which ions will start being "too light" to represent pro- 
tons? 

Even before answering this question, one could ask 
whether the PIC simulation plasmas could be governed 
by some similarity laws involving the mass ratio. Similar- 
ity theory has been applied successfully in hydrodynam- 
ics. It allows us to predict certain properties of an object 
from experiments performed with its miniaturised model. 
Well-known cases of such experiments involve pumps, 
turbines or aircrafts [37], [HI . Similarity laws have also 
been derived for magnetic confinement fusion (see [39| , 
and references therein) or relativistic laser-plasma inter- 
actions and laboratory astrophysics |40|-|42|. Similarity 
laws would allow us to compensate a reduced mass ratio 
with some other parameter, by which the computational 
efficiency can be altered. 

We show in section 2 that it is not possible to derive an 
universal description of the growth rate spectrum, which 
is not explicitly dependent on the mass ratio. Section 3 
will therefore aim at providing a restricted solution to the 
problem. PIC simulations typically probe the long term 
nonlinear evolution of an unstable beam-plasma system. 
The unstable spectrum usually contains more that one 
unstable mode, and these modes grow at different ex- 
ponential rates. A hierarchy of unstable modes can be 
established in terms of their growth rate, and a crite- 



rion can be imposed on the mass ratio by demanding 
that this hierarchy be preserved. The consequences of 
this condition are then explicitly calculated for the case 
of a cold relativistic electron beam passing through a 
un- magnetized and cold plasma [43|. Section 4 is the 
discussion, which brings forward a possible explanation 
why the shock formation in Ref. Q does apparently not 
depend on the mass ratio. 

II. ABSENCE OF A SIMILARITY LAW 

Consider a problem that involves n independent di- 
mensional variables (xi, . . . , x n ), and m fundamental di- 
mensions such as meter, second, etc. The so-called Buck- 
ingham's method |44[ to reduce the number of variables 
and to derive similarity laws is the following. 

Identify the pairs of x^s that share the same physi- 
cal unit. If this is the case for the variables x kl and 
x k2 , then replace (x kl ,x k2 ) by (x kl ,x kl /x k2 ) in the list 
of variables. After iterating this process for any such 
pairs, we are left with the modified set of variables 
Oi, . . . /xk 2 , • • • ,x kl _ 1 /x kl ) where I + m = n. 

Buckingham's "II theorem" then states that any un- 
known function of the form 

ffa, . . . ,x m ,x kl /x k2 , . . . ,x kl _Jx kl ) = 0, (1) 

can indeed be cast under the form, 

0(7Ti, . . .,-Km- p ,X kl /x k2 ,. . .,X kl _Jx kl ) = 0, (2) 

where the variables i\i are dimensionless products of the 
initial a^'s, and p the number of fundamental dimensions 
among (xi, . . .,x m ). 

Consider as a simple illustration a swinging pendu- 
lum with the mass M (kg) and the length I (m), which 
oscillates with the constant period T (s) in a gravita- 
tional field g (m/s 2 ). The first step of Buckingham's 
method, namely the pairing of variables that share the 
same dimension, can be skipped here since all 4 variables 
(M, Z,T, g) have different dimensions (kg, m, s, m/s 2 ). 
We thus have here m = 4 (4 variables) and only p = 3 
(kg, m, s) as g does not add any extra fundamental di- 
mension to the problem. 

Buckingham's theorem states in this case that 
any function f(M,l,T : g) = can be expressed as 
(/>(7T m _ p =i) = 0, so that the problem is eventually a func- 
tion of one single dimensionless parameter. A subsequent 
dimensional analysis shows that the universal parameter 
must be a power of gT 2 /L. Clearly, the mass M cannot 
participate in the dimensionless parameter, because no 
other variables could cancel its physical unit. The period 
T is thus independent of the mass M and only a func- 
tion of the ratio g/L. Note that the theorem does not 
distinguish between "input" variables (what is known, 
e.g. M, l,g) and "output" variables (what is looked for, 
e.g. T). Each quantity is treated in the same way and 
all contribute to m and to p. 
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Turning now to the present problem, we see from the 
first step of Buckingham's method that a similarity law 
without an explicit dependence on the mass ratio can- 
not exist. Whatever the list of variables describing the 
problem may be, the mass parameters (m e ,M p ) will be 
a part of it. The first step of the process will just re- 
place (m e ,Mp) by (m e ,m e /M p ), and Buckingham's the- 
orem reduces the number of variables left once all the 
dimensionless trivial ratios have been formed. At any 
rate, Buckingham's theorem states that the mass ratio 
remains as an explicit parameter in the final reduced set 
of parameters. 

There is therefore no hope of unraveling similarity laws 
connecting two systems (A) and (B) with different mass 
ratios. A change of the mass ratio cannot be "compen- 
sated" by a shift of the other variables. Buckingham's 
analysis of the problem proves an intuitively simple rea- 
soning: electrons and ions define different time scales in 
terms of their respective mass. The time evolution of 
the system can be normalized to any one of them, but 
it cannot fit both at the same time. Note that although 
the rest of the article focuses on the linear regime of an 
electron beam plasma system, the present conclusion is 
very general, and valid for the overall evolution of any 
kind of system comprised of two species. 

Let us initially assume that the ions are immobile, 
yielding an electron-to-ion mass ratio of zero. Electrons 
are therefore the only population bringing a mass into the 
parameter list. The electron mass m e will therefore ap- 
pear among the x ]_ ... x m Eq. (pQ). Buckingham's the- 
orem here states that these m dimensional variables can 
be replaced by m— p dimensionless variables 7r± . . . 7r m _ p . 
Because no mass ratio can appear among the I — 1 ratios 
which are arguments of the function <\> in Eq. (j2j) , the un- 
derlying equations can not rely explicitly on the electron 
mass. Equation (j2j) shows that m e must have been "ab- 
sorbed" by one of the dimensionless 7r's. This is precisely 
what is observed when dealing with such questions: the 
time parameter is frequently normalized to the electronic 
plasma frequency which includes the electron mass. 

III. PRESERVING THE GROWTH RATE 
HIERARCHY OF THE UNSTABLE MODES 

In ultra-relativistic laser-plasma interaction, similarity 
theory states that laser-plasma interactions with different 
ao = eAo/m e c 2 and n e /n c are similar, as soon as the 
similarity parameter S = n e /aon c is the same [42[ (Aq 
is the laser amplitude, n e is the plasma electron density, 
and n c = ra e o;o/47re 2 is the critical density for a laser 
with frequency c^o- 

The previous Buckingham analysis has demonstrated 
that there cannot be any such similarity parameter in 
the problem we consider here. As soon as the ions are 
allowed to move, the electron to ion mass ratio R must 
appears explicitly in any list of dimensionless parameters 
describing the system. Two systems differing only by 



their mass ratio will not evolve similarly. 

A deviation of the simulation dynamics from the true 
plasma dynamics is acceptable, as long the modifications 
are only quantitative. The simulation can in this case 
still provide important qualitative insight into the plasma 
evolution, which can not be obtained by any other means. 
However, somewhere in between the mass ratio of 1/1836 
and the (positron) mass ratio of 1, a line must be crossed 
when even this is not the case any more. 

For an unstable system with ratio R = ra$/ m e between 
the ion and the electron mass, the unstable spectrum 
§ = {k £ M 3 I <S(k, R) > 0} is comprised of all the modes 
with wavenumbers k with an amplitude that grows at 
the exponential (positive) growth rate S. Among these 
modes, the most unstable mode k m (i2) defined by, 

<S(k m , R) = max{£(k, #)} k es S m (R), (3) 

plays a peculiar role because it is the one which growth 
determines the outcome of the linear phase. The evolu- 
tion of k m , as a function of R, may be continuous, or 
not. For clarity, let us consider the example studied in 
Sec. IIII A[ of a one dimensional beam-plasma system. 
The dispersion equation in this case gives two kinds of 
unstable modes: the two-stream modes and the Buneman 
modes. Let us assume that we have plasma parameters 
that are such that the dominant mode k m is a two-stream 
mode. It is possible to find a range of mass ratios R, for 
which the two-stream mode always grows fastest. In this 
case, k m evolves continuously with R. A variation of the 
mass ratio may trigger in another case a transition from 
the two-stream regime to a Buneman regime where k m 
is the wavenumber of a Buneman mode. Here, the evolu- 
tion of k m will be discontinuous and we will talk about 
an altered mode hierarchy. 

Changing the mass ratio in such a way that the mode 
hierarchy is altered will thus result in a different plasma 
evolution during the linear growth phase of the instabil- 
ities. In our ID example, the typical size of the patterns 
generated in the early evolution will change abruptly by 
a factor n&/n e , where n^e are the beam and plasma elec- 
tronic densities respectively. For the 2D system consid- 
ered in Sec. IIII Bl a switch from a two-stream to a fila- 
mentation regime results in the generation of magnetized 
filaments instead of electrostatic stripes. 

Although we focus in what follows on a specific setup, 
most kinds of beam-plasma systems encountered in the 
literature also exhibit more than one type of unstable 
modes @,|H,|46|. 

The criterion we propose is that the modified mass ra- 
tio must not alter the mode hierarchy. Note that this is a 
necessary but not a sufficient condition. If the mode hier- 
archy is altered, then the evolution of the system should 
be affected as well. However, even a similar linear phase 
could result in a different non-linear long-term evolution 
prompted by a different mass ratio. 

Even if thermal effects are neglected, the analysis of 
the full spectrum of unstable waves is involved for ener- 
getic astrophysical plasmas [46| . The identification of the 



fastest-growing wave mode requires the evaluation of the 
full three-dimensional spectrum of wave vectors 
I49I ]. In what follows, the proposed criterion for the mass 
ratio is applied to the simple and generic, yet important, 
system formed by a relativistic electron beam that passes 
through a plasma with an electronic return current. The 
return current initially cancels the beam current and the 
ion charge density cancels the total electronic one. Since 
some PIC simulations are still performed in a ID geome- 
try, we start by analyzing the ID case before we turn to 
the more realistic 2D and 3D ones. 



A. Relativistic electron beam - ID simulation 

We consider a relativistic electron beam with the den- 
sity rib •> the velocity V5 and the Lorentz factor 75 = 
(1 — v 2 /c 2 ) -1 / 2 , which passes through a plasma with 
the ion density rii and the electron density n e with 
Tii = rib + n e . The drift velocity v e of the electrons of 
the background plasma fulfills ribVb + ^ e v e = 0. For the 
stability analysis, we consider the response of the system 
to harmonic perturbations ex exp(zk-r — iuot). We assume 
that the particles move along the z-axis and we consider 
only wavevectors with k || z. All plasma species are cold. 
The dispersion equation is readily expressed as the sum 
of the contributions by the beam, by the return current 
and by the ions with mass Mi. 
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where 7 e = (1 — v 2 /c 2 ) -1 / 2 is the Lorentz factor of the 
return current. We introduce with uo 2 = Aim e e 2 /m e the 
dimensionless variables 
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The dispersion equation expressed in these variables is 
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As long as a <C 1 the return current remains non- 
relativistic with j e ~ 1. For the strictly symmetric case 
with a = 1, the return current becomes relativistic with 
7 e = 76 . The dispersion equation (|6)) defines two kinds of 
unstable modes [50| . The two-stream instability is driven 
by the two electron beams. In the limit a <C 1 this insta- 
bility has its maximum growth rate 5 at the wavenumber 
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FIG. 1. (Color online) Growth rates of the two branches of 
unstable electrostatic modes derived from Eq. © for a = 0.2, 
R = 1/1836 and 75 = 4. Two-stream modes reach their 
maximum growth rate for Z ~ 1 and Buneman modes for 
Z - 1/a. 
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FIG. 2. Separatrices of the parameter space intervals domi- 
nated either by the two-stream instability or by the Buneman 
instability. The curves R m correspond to different mass ra- 
tios for the electrostatic ID system considered in Section lHI Al 
Plain lines: numerical evaluation. Dashed lines: borders de- 
fined by Eq. ©. 



Figure [T] displays the growth rate curves obtained from 
Eq. ([6j). Both wave branches share the same Z-interval 
if a ~ 1. Equations (|7I8|) show how the mode hierarchy 
relies explicitly on the mass ratio. Only the growth rate 
of the Buneman modes scales like R 1 / 3 . Changing the 
mass ratio can thus change the mode hierarchy. 

Figure [2] depicts the range of parameters (75, a) where 
two-stream and Buneman modes govern the spectrum for 
various mass ratios R. The separatrix between the do- 
main is plotted for R = 1/30, 1/100, 1/400 and 1/1836. 
The Buneman modes grow faster below the curve, while 
the two- stream modes are dominant above. The separa- 
trix Rm between the two domains is given by 



The unstable Buneman modes arise from the interaction 
of the electronic return current with the ions. These ad- 
ditional modes grow for a <C 1 at the wavenumber 
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Z ~ 1/a, with 5 ~ -^i? 1/3 . 
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Let us assume that we run a ID PIC simulation from 
the parameters pictured by the circle labeled "PIC1". 
For R = 1/1836, the corresponding system lies in the 
two-stream region. As we increase R, the growth rate 
of the Buneman instability increases relative to that of 
the two-stream instability. For sufficiently light ions, the 
Buneman instability can even outgrow the two-stream 
instability. Given some simulation parameters (75, a), 
the largest mass ratio that leaves the mode hierarchy un- 
changed is readily calculated from Eq. (|9]) if a <C 1. One 
can see that a mass ratio as high as 1/30 is allowed only 
for weakly relativist ic systems. If a PIC simulation uses 
the parameter values denoted by "PIC2" , the present cri- 
terion does not restrict the mass ratio. Any value larger 
than 1/1836 would be in favor of the Buneman modes, 
which are already governing the system. 

B. Relativistic electron beam - 2D and 3D 
simulations 




FIG. 3. (Color online) Growth rate map as a function of 
Z x , Z z . The parameters are the same as in Fig. [2] The beam 
velocity vector points along the z axis. 



The previous reasoning is now expanded to a 2D geom- 
etry. It is equivalent to a full 3D geometry with regard 
to a linearized theory, as long as the system is cold and 
does not have two distinct symmetry axes. An exam- 
ple is a beam velocity vector that is not aligned with 
the magnetic field direction. Here, the v& forms the sole 
symmetry axis and it defines one direction. A second 
dimension takes into account the unstable modes with 
wave vectors that are not parallel to V5. These modes 
compete with the two-stream and Buneman modes. One 
finds the Weibel (or filamentation) modes for k _L v&, 
which could play a major role in the magnetic field gen- 
eration that is necessary to explain Gamma Ray Bursts 
[H, HO). For obliquely oriented wave vectors, the so-called 
"oblique modes" are likely to govern parts of the rela- 
tivistic regime p52j . 

The dispersion equation is more involved in 2D than in 
ID, because unstable modes are generally not longitudi- 
nal (i.e. electrostatic with k || E). While oblique unsta- 
ble modes have been known to exist for some decades now 
[53l-[55|, the first exact cold fluid analysis of the full un- 
stable spectrum was only recently performed by Califano 
et. al. (43|| . The dielectric tensor is computed exactly 
from the Maxwell's equations, the continuity equation 
and the Euler equation for the three species involved. 
We choose v& || z and a wave vector (k x ,0,k z ) in the 
(x,z) plane. The normalized wave vector Z from Eq. ([5]) 
is now extended to two dimensions, 



The dielectric tensor has been computed symbolically us- 
ing a Mathematica Notebook designed for this purpose 
[56j. The dispersion equation reads now 

det(T) = 0, (11) 



where the tensor T is specified in the Appendix. The 
dispersion equation is an 8th degree polynomial. Its nu- 
merical solution is straightforward. Figure [3] shows a plot 
of the growth rate in terms of Z = (Z x , Z z ) for the same 
parameters as Fig. [TJ The Weibel or filamentation modes 
are characterized by Z z = 0, the two-stream and Bune- 
man modes by Z x =0, while the oblique modes consti- 
tute the remaining spectrum. The ridge in the growth 
rate map at low Z z stems from the interaction of the two 
electron beams, while the interaction of the ions with the 
bulk electrons is responsible for that at larger Z z s. 



The growth rates of the Weibel modes and of the 
oblique modes can be estimated for immobile ions with 
R = and for low a as [55| , 

Weibel: S w = $J^- 
V 76 

Oblique: *o = ^{^'\ (12) 

These expressions must be corrected in a nontrivial way 
in the ultra-relativistic limit as a approaches unity. Prior 
to the formulation of our criterion for the mass ratio 
i?, we elucidate the hierarchy map and how it evolves 
with R. Figure |4] pictures the separatrices of the do- 
mains in the parameter space for R = 1/1836 and for 
R = 1/30. Weibel modes tend to govern the high 
beam density regime, and slightly expand the mildly- 
relativistic (around 75 = 20) part of their domain when 
R grows. Buneman modes modes govern the lower-right 
corner of the graph, and being scaled like i? 1 / 3 , increase 
their domain as well with R. As a result, the domains 
governed by the oblique modes shrinks with a growing R. 
The two-stream modes actually never govern the system, 
because they are outgrown by the oblique modes as soon 
as 76 > 1. 
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FIG. 4. (Color online) 2D hierarchy map in terms of (75, a). 
Plain lines: R = 1/1836. Dashed lines: R = 1/30. Weibel 
instability tend to govern the high density regime, Buneman 
the ultra-relativistic one, and oblique the rest of the phase 
space. 

The parameter space diagram reveals a triple point, at 
which the separatrices merge. For R = 1/30, its coor- 
dinates are (75, a) ~ (30,0.48). Figure [5] shows how the 
triple point location evolves towards the ultra-relativistic 
regime 75 ^> 10 3 as the ion mass is increased to that of 
a proton. This triple point is not likely to be important 
in astrophysical flows, since even the Lorentz factors of 
GRB jets do not reach such high values. However, it can 
become an issue in PIC simulations, where mass ratios of 
30 and Lorentz factors of a few tens are not uncommon. 
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FIG. 5. (Color online) Coordinates of the triple point where 
oblique, Weibel and Buneman modes grows exactly the same 
rate, in terms of the electron to proton mass ratio. 
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We now compute the largest R that leaves unchanged 
the mode hierarchy for a given parameter set (75, a), and 
display the result on Fig. [6] (in fact, the smallest in- 
verse mass ratio i?" 1 is plotted, for better clarity). If, 
for R = 1/1836, a system lies in the Buneman region, 
then increasing R will not change the dominant mode. 
The same holds for the systems pertaining initially to 
the Weibel zone. But the system represented by "PIC1" 
on Fig. [4] remains governed by the oblique modes only 
up to a certain value of the mass ratio, beyond which it 
goes over into the Buneman domain. The same can be 
said for "PIC2" : initially lying in the oblique domain, it 
goes over into the Weibel domain beyond a critical value 
of the mass ratio. Only systems already located in the 
oblique domain for R = 1/30 continue to do so as we alter 
R from 1/1836 to 1/30. Of course, we speak here only 



FIG. 6. (Color online) Smallest inverse mass ratio value R^ 1 
leaving unchanged the 2D modes hierarchy for a given param- 
eter set (75, a) for 1/1836 < R m < 1/30. The uniform white 
region refers to configuration where the present criterion does 
not constrain the mass ratio. 

about the lower part of the graph that corresponds to 
small values of a. The dominant mode depends through 
Rmijb,®) critically and in a nontrivial way on both, a 
and on 75 in the upper part of Fig. [6j 

We thus find a significantly extended region of the pa- 
rameter space, namely, the uniform white domain on Fig. 
[6j where the criterion that the mode hierarchy be un- 
changed does not restrict the value of the mass ratio. For 
a system lying in this region, the dominant mode is the 
same, regardless of whether R = 1/1836 or 1/30. In the 
lower-right corner (i.e., diluted ultra-relativistic beams), 
the border is defined by the equality of the oblique mode 
(see Eq. [12]) with the Buneman one for R = 1/1836, 

<* = #76|i?=l/1836- (13) 
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In the lower-left corner (i.e, diluted, weakly relativistic 
regime), the border is determined by equating the oblique 
growth rate with the Buneman one, but now for R = 
1/30, 

OL = #76|i?=l/30. (14) 

The upper-border of the uniform Weibel domain is an- 
alytically more intricate. Let us just mention that the 
particular shape exhibited for 75 ~ 2 arises from the 
Weibel growth rate which reaches a maximum around 
this value. Expression ([T2]) for this quantity makes it 
clear that 5w(vb — 0) = 0, while lim v? ^ c 8w = 0. As a 
consequence, Sw reaches a maximum for an intermediate 
Lorentz factor 75 = a/3, which is easily calculated from 
Eq. ([T2]) . Although this value is not exact for a close to 
unity and for R ^ 0, the "Weibel optimum" for 75 stays 
close to a/3, explaining the "bump" at this location. 

IV. CONCLUSION 

The importance of the electron to ion mass ratio R = 
m e /mi for the realism of PIC simulations has been ad- 
dressed here from an analytical point of view. Since there 
cannot be any rigorous similarity theory encompassing 
this quantity, an attempt has been made to identify a 
threshold R mi beyond which a given simulation can no 
longer be trusted to be physically accurate during the 
initial exponential growth phase of the instability. This 
initial wave growth can be addressed by a linearized the- 
ory. 

Whether it be relevant for astrophysical plasmas or for 
inertial fusion, many systems investigated through PIC 
simulations give rise to the growth of waves that can be 
addressed by an analysis of the linear dispersion relation. 
The idea is therefore to find the maximum value of R, 
which leaves unchanged the hierarchy of the linearly un- 
stable modes. The condition we propose is necessary but 
not sufficient: for the system evolution to be preserved, 
the linear evolution and, more specifically, the type of 
the fastest growing mode must remain unchanged as we 
change R. However, two similar linear growth phases can 
eventually result in a different non-linear state. 

Because the application of the criterion depends on the 
linear unstable spectrum, and therefore on the system 
under scrutiny, we have focused on the generic system 
formed by a relativistic electron beam passing through 
a plasma with return current. For a ID simulation, 
the competing modes are the two-stream mode and the 
Buneman mode (see Fig. [2]). The criterion of the pre- 
served mode hierarchy does provide a value of R mi if the 
spectrum is governed by the two-stream instability for 
R = 1/1836. As a result, for example, the simulation 
of a 10 times diluted beam with 75 = 2 cannot be per- 
formed with ions that have a mass below ~ 100 times 
heavier than that of the electrons. For systems governed 
by the Buneman instability when R = 1/1836, our crite- 
rion does not give any upper value of R m . 



The 2/3D case is even more interesting as more modes 
compete in the linear phase. Here, the Buneman, the 
oblique and the Weibel instabilities can dominate the 
linear phase, while the two-stream instability is unim- 
portant for relativistic beam speeds. For R = 1/30 and 
1/1836, the hierarchy map is plotted on Fig. [3]in terms of 
the density ratio a and the beam Lorentz factor 75. One 
can notice how the upper part (a ~ 1) does not vary with 
R. This could explain why PIC simulations of collisions 
between equally dense plasma shells did not show much 
difference as the mass ratio has been altered [4]. The 
dominant mode is definitely the Weibel (filamentation) 
one in this region. Things should be different when sim- 
ulating collisions of shells with a different density, like in 
[I?], [Hj, and varying the mass ratio. 

The value of R m in terms of (75, a) is predictable at 
low a, and more involved for a ~ 1 (see Fig. [6]). A 
few points can be emphasized at this junction: First, the 
most sensitive points are the ones located near a border 
between two modes for R = 1/1836. When R departs 
from this value, the border moves, say from mode (A) 
domain to mode (B) domain, so that (A) increases to 
the expenses of (B). If the point was initially in the (A) 
domain, it remains there and the criterion is not bind- 
ing. But if the point was close to the border, yet in the 
(B) region, then a slight increase of R transfers it to the 
(A) region. This is why on Fig. [6j the white region of 
unconstrained R always borders Rm = 1/1836. 

Second, some alteration of the mode hierarchy are 
more dramatic than others. Figures |4] & [6] show that 
3 kinds of transitions can be triggered when increas- 
ing R: oblique to Buneman (OB - for diluted beam), 
oblique to Weibel (OW - high density, weakly rela- 
tivistic) and Buneman to Weibel (BW - high density, 
ultra-relativistic). For diluted beams, the OB transi- 
tion switches the wavelength of the dominant mode from 
Z z ~ 1 to 1/a, resulting in generated structures a times 
smaller. Furthermore, oblique modes generate partially 
electromagnetic transverse structures whereas the elec- 
trostatic Buneman modes do not. More dramatic can 
be the OW transition as we now switch from a quasi- 
electrostatic dominant modes to an electromagnetic one. 
But the BW transition is by far the most powerful as 
the generated patterns switch from stripes (Buneman) 
to Filaments (Weibel). 

Note however that transitions are smoother than they 
appear because the switch from one mode regime to an- 
other is not immediate when a border is crossed. Suppose 
we move from domain (A) to (B). As we approach the 
border, mode (A) keeps growing faster, but mode (B) 
grows almost as fast, until the growth rates are strictly 
equal right on the border. There is therefore a zone ex- 
tending on both side on the line where a proper inter- 
pretation of the linear regime needs to account for the 
growth of (A) plus (B), thus smoothing out the transi- 
tion. 

The present article is a first step towards a system- 
atic search. The method proposed has been applied to a 
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generic beam-plasma system, evidencing non-trivial val- 
ues of R m . A similar analysis can be easily conducted 
varying the set-up: one needs first to evaluate the growth- 
rate map (the counterpart of Fig. [3j) as a function of k 
for the system under scrutiny with R = 1/1836. The 
same plot is then evaluated for the desired value of R. If 
the dominant mode remains the same, then the present 
criterion is met. 

For example, when dealing with the problem of mag- 
netic field amplification and particles acceleration in Su- 
pernova Remnants, a typical PIC setup consists in a non- 
relativistic beam of protons passing through a plasma 
with a guiding magnetic field [HI, [59| . Due to the mag- 
netization, unstable modes such as the Bell's ones [45j 
enrich the spectrum, and it would be interesting to ap- 
ply our criterion also to these cases. 
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Appendix A: 2D and 3D Tensor 

Choosing the axis z for the flow direction, and Z = 
(Z x , 0, Z z ), the tensor involved in the dispersion equation 
(fTTj) is symmetric, and reads, 



T 



Tn T 3 i 
T 22 
Tsi T33 



(Al) 



where, 
Tu = l 

T22 = — 
T 33 = l 



R(l + a) 1 

x [ 



a 1 

lb 7e 



2 

+ Z 2 ) 7 6 + <*7e + (R - x 2 + Ra) 75 



x 2 (3 2 
R(l + a) 



X 2 Jble 



1 



(x-Z z ) 2 jI (x + Zza) 2 ^ 

„2 



or 



x 2 [f3 2 (x-Z z ) 2 lb (x + Z z a) 2 j e 



T31 — — T 



P 2 XJ e + Z z aj e [Z z - X)j b 
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